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Abstract. We present a continuous finite element method for some examples 
of fully nonlinear elliptic equation. A key tool is the discretisation proposed in 
Lakkis & Pryer (2011) allowing us to work directly on the strong form of a lin- 
ear PDE. An added benefit to making use of this discretisation method is that 
a recovered (finite element) Hessian is a biproduct of the solution process. We 
build on the linear basis and ultimately construct two different methodologies 
for the solution of second order fully nonlinear PDEs. Benchmark numeri- 
cal results illustrate the convergence properties of the scheme for some test 
problems as well as the Monge— Ampere equation and the Pucci equation. 
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1. Introduction 

Fully nonlinear PDEs arise in many areas, including differential geometry (the 
Monge- Ampere equation), mass transportation (the Monge-Kantorovich problem), 
dynamic programming (the Bellman equation) and fluid dynamics (the geostrophic 
equations) . The computer approximation of the solutions of such equations is thus 
an important scientific task. There are at least three main difficulties apparent to 
someone attempting to derive numerical methods for fully nonlinear equations: first, 
the strong nonlinearity on the highest order derivative which generally precludes a 
variational formulation, second, a fully nonlinear equation does not always admit 
a classical solution, even if the problem data is smooth, and the solution has to 
sought in a generalised sense (e.g., viscosity solutions), which is bound to slow 
down convergence rates, and third, a common problem in nonlinear solvers, the 
exact solution may not be unique and constraints, such as convexity requirements 
must be included in the constraints to ensure uniqueness. 

Regardless of the problems, the numerical approximation of fully nonlinear sec- 
ond order elliptic equations^ as described in Caffarelli and Cabre 1995 , have been 
the object of considerable recent research, particularly for the case of Monge- 
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For more general classes of fully nonlinear equations some methods have been 
presented, most notably, at least from a theoretical view point, in |B6hmer| |2008 



where the author presents a finite element method shows stability and consis 
tency (hence convergence) of the scheme, following a classical "finite difference" 
approach outlined by Stetter] 1973 which requires a high degree of smoothness on 
the exact solution. From a practical point of view this approach presents difficulties, 
in that the finite elements are hard to design and complicated to implement, 
in Davydov and Saeed 120101 a useful overview of Bezier-Bernestein splines in two 



spatial dimensions is provided and a full implementation in [Davydov and Saeed| 
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2012 . Similar difficulties are encountered in finite difference methods and the con- 



cept of wide-stencil appears to be useful, for example by Kuo and Trudinger 1992 
20^5) |Oberman| [2008' , F roesel 
In |Feng and Neilan, [2009a|b 



2011 



Awanou 2010 the authors give a method in which 



they approximate the general second order fully nonlinear PDE by a sequence of 
fourth order quasilinear PDEs. These are quasilinear biharmonic equations which 
are discretised via mixed finite elements, or using high-regularity elements such as 
splines. In fact for the Monge-Ampere equation, which admits two solutions, of 
which one is convex and another concave, this method allows for the approximation 
of both solutions via the correct choice of a parameter. On the other hand although 
computationally less expensive than finite elements (an alternative to mixed 
methods for solving the biharmonic problem) , the mixed formulation still results in 
an extremely large algebraic system and the lack of maximum principle for general 
fourth order equations makes it hard to apply vanishing viscosity arguments to 
prove convergence. A somewhat different approach, based on C'^-penalty, has been 



2011 



recently proposed by Brenner et al. 2011a , as well as "pseudo time" one by Awanou 



It is worth citing also a least square approach described by |Dean and Glowinski| 
This method consists in minimising the mean-square of the residual, using 



2006 



a Lagrange multiplier method. Also here a fourth order elliptic term appears in the 
energy. 

In this paper, we depart from the above proposed methods and explore a more 
"direct" approach by applying the nonvariational finite element method, introduced 

as a solver for the Newton iteration directly derived 



m 



Lakkis and Pryer 2011 



from the PDE. To be more specific, consider the following model problem 

^[u] := F{D'^u) - / = 



(1.1) 



with homogeneous Dirichlet boundary conditions where / : i? — > IR is prescribed 
function and F : Sym(IR'^^'*) — > R is a real- valued algebraic function of symmetric 
matrixes, which provides an elliptic operator in the sense of |Caffarelli and Cabre| 
1995 , as explained below in Definition 2.2 The method we propose, consists in 



applying a Newton's method, given below by equation (4.2 ) of the PDE ( |1.1[ ), which 
results in a sequence of linear nonvariational elliptic PDEs that fall the framework 



of the nonvariational finite element method (NVFEM) proposed in ILakkis and 



Pryer 2011 . The results in this paper are computational, so despite not having 
a complete proof of convergence, we test our algorithm various problems that are 
specifically constructed to be well posed. In particular, we test our method on the 
Monge-Ampere problem, which is the de-facto benchmark for numerical methods 
of fully nonlinear elliptic equations. This is in spite of Monge-Ampere having an 
extra complication, which is conditional ellipiticity (the operator is elliptic only 
if the function is convex or concave. A crucial, empirically observed feature of 
our method is that the convexity (or concavity) is automatically preserved if one 
uses elements or higher. For elements this is not true and the scheme must 
be stabilized by reenforcing convexity (or concavity) at each timestep. This was 
achieved in Pryer 2010 using a semidefinite programming method. In a different 
spirit, but somewhat reminiscent, a stabilization procedure was obtained in |Breniier] 



et al. 2011a by adding a penalty term. 



The rest of this paper is set out as follows. In S|2]we introduce some notation, 
the model problem, discuss its ellipticity and Newton's method, which yields a 
sequences of nonvariational linear ised PDE's. In ^ we review of the nonvariational 
finite element method proposed in|Lakkis and Pryer 2011 and apply it to discretise 
the nonvariational linearised PDE's in Newton's method. In ^ we numerically 
demonstrate the performance of our discretisation on a class of fully nonlinear 
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PDE, those that are elhptic and well posed without constraining our solution to 
a certain class of functions. In fjS] we turn to conditionally elliptic problems by 
dealing with the prime example of such problems, i.e., Monge-Ampere . We apply 
the discretisation to the Monge-Ampere equation making use of the work [Xguilera| 
and MorinI |2009| to check finite element convexity is preserved at each iteration. 



Finally in ^|6] we address the approximation of Pucci's equation, which is another 
important example of fully nonlinear elliptic equation. 

All the numerical experiments for this research, were carried out using the 
DOLFIN interface for FEniCS Logg and Wells 2010 and making use of Gnuplot 
and Para View for the graphics. 



2. Notation 



2.1. Functional set-up. Let J7 C R'' be an open and bounded Lipschitz domain. 
We denote L2(/2) to be the space of square (Lebesgue) integrable functions on 
f2 together with its inner product {v,w) := J^vw and norm \\v\\ ll^llL2(f2) ^ 

1/2 

{v,v) . We denote by {v \ w) the action of a distribution v on the function w. 

We use the convention that the derivative Dm of a function m : ]7 — >■ R is a row 
vector, while the gradient of u, Vm is the derivatives transpose (an element of R'^, 
representing Dm in the canonical basis). Hence 



Vm = (Dm)"^. 



(2.1) 



For second derivatives, we follow the common innocuous abuse of notation whereby 
the Hessian of u is denoted as D^u (instead of the more consistent DVm) and is 
represented by a ci x d matrix. 

The standard Sobolev spaces are Ciarlet [1978 , Evans 1998 



H'=(I?) := W^(I?) 



e L2(I2) : D°0e L2(r2) \ , (2.2) 

|a|<fc J 

Rlin) closure of C^{f2) in (2.3) 



where a = {ai, ad} is a multi-index, \a\ = '^j ^^"^ derivatives D" are 

understood in a weak sense. 



We consider the case when the model problem (1.1 1 is uniformly elliptic in the 
following sense. 



2.2. Definition (ellipticity |Caffarelli and Cabri] [1995| ). The operator ^[ J 
in Problem ( 1.1 ) is called elliptic on'^ C Sym (R"^"*) if and only if for each M e '^S' 
there exist A > A > 0, that may depend on M such that 

A sup |iV|| < i^(Af + AT) - F(M) < ^ sup lAT^I V TV e Sym (R'^^''). (2.4) 



If the largest possible set ^la for which (2.4) is satisfied is a proper subset of 
Sym(R''^'^) we say that ^[-J is conditionally elliptic. 



The operator .yV[-] in Problem ( 1.1 1 is called to be uniformly elliptic if and only 



if for some A, yl > 0, called ellipticity constants, we have 



A sup |Af^| <F(M + Ar)-F(M) <yl sup |7V^| V7V, M e Sym (R"^'*). (2.5) 
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If F is differentiablc (2.5) can be obtained from conditions on the derivative of 



F. A generic M e Sym (R''^'') is written as 



M 



jrid,! ■ ■ ■ md,d_ 
so the derivative of F in the direction TV is given by 

BF{M)N = F'{M):N 
where the derivative matrix F'{M) is defined by 

'dF{M) /dm^ ^ ... dF{M) /dm^ , 



F'{M) := 



dF{M)/dm^^^ 



dF{M)/dm,, 



(2.6) 



(2.7) 



(2.8) 



Suppose F is differentiablc. Then (2.4) is satisfied if and only if for each M € 
there exists > such that 

F' {M)i> V^eR''. (2.9) 

Furthermore = Sym(IR''^'') and /i is independent of M if and only if (2.5) is 
satisfied. 

2.3. Assumption (smooth elliptic operator). In the remainder of this paper 
we shall assume that J^W is conditionally elliptic on '€ and 



F e c^C^). 

Unless otherwise stated we will also assume that 



(2.10) 



Sym ( 



Ddxd\ 



2.4. Newton's method. The smoothness assumption 2.3 allows to apply New- 
ton's method to solve Problem (1.1). 

Given the initial guess u'^ € C (i7), with D^ii° e for each n e Nq, find 

^n+l g (j2^^^ ^j^j^ D2m"+1 e such that 

D^[m"] - u") = -^K], (2.11) 

where Dc/K[u] indicates the (Frechet) derivative, which is formally given by 

.yV[u + €v] - jV[u] 



B.yi^[u]v — lim 



e->0 e 

F(D^u + eD^v)- FiB^u) 
~ Imi 



(2.12) 



= F'(D^u) : Wv, 



for each v G C^(/2). Combining (2.11) and (2.12) then results in the following 
nonvariational sequence of linear PDEs. Given u^for each rt S Nq find m""^^ such 
that 

(2.13) 



71 + 1 



/-F(DV") 



The PDF (2.13) comes naturally in a nonvariational form. If we attempted 
to rewrite into a variational form, in order, say, to apply a "standard" Galerkin 
method, we would introduce an advection term which would depend on derivatives 
of F' , i.e., for generic v, w 



F'{D^v):J:)^w = div[F'{D^v)Vw] - div[F'{D^v)] Vi 



(2.14) 
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where the matrix-divergence is taken row-wise: 



and the chain rule provides us, for each j = I, . . . , d, with 



(2.15) 



d 



(2.16) 



k,l=l 



This procedure is undesirable for many reasons. Firstly it requires F to be twice 
differentiable and it involves a third order derivative of the functions and 
m" appearing in (2.111. Moreover, the "variational" reformulation could very well 
result in the problem becoming advection dominated and unstable for conforming 



FEM, as was manifested in numerical examples for the linear equation Lakkis and 



Pryer 2011 §4.2]. In order to avoid these problems, we here propose the use of the 



nonvariational finite element method described next. 



3. The nonvariational finite element method 



In Lakkis and Pryer 2011 we proposed the nonvariational finite element method 



(NVFEM) to approximate the solution of problems of the form (2.13). We review 



here the NVFEM and explain how to use it in combination with the Newton method 
to derive a practical Galerkin method for the numerical approximation of Problem 



( 1.1 1's solution. 



3.1. Distributional form of (2.13) and generalised Hessian. Let A £ Loo(^)'^^'' 
and for each x G f2, let A{x) G Sym (R''^'^), the space of bounded, symmetric, pos- 
itive definite, d x d matrices and / : i7 — > R. The Dirichlet linear nonvariational 
elliptic problem associated with A and / is 

A: D^u = / and = 0. (3.1) 

Testing this equation, and assuming u G H^(J7)nlIo(/2) such that Vu|g^ G L2{df^), 
we may write it as 



{A:D^u,cf>) = {f, 



(3.2) 



To allow a Galerkin type discretisation of (3.2 1, we need to restrict the test functions 
(/) to finite element function spaces that are generally not subspaces of II^(i7). So 
before restricting, we need to extend and we use a traditional distribution-theory (or 
generalised-functions) approach. Given a function v G Y{^{fl) and let n : dQ — > 
be the outward pointing normal of fi then the Hessian of w, TP'v satisfies the 
following identity: 



(D^t;,^) = - / Vf ® V(/)+ / \/v®n(j) V(/)GH1(/?), 
Jn Jdn 



(3.3) 



where a®b:= ab^ for a, 6 column vectors in R"^. If u G II^(!7) with Vti|g^ G 



H ^^"^{dn) the right-hand side of (3.3) still makes sense and defines D^w as an 
element in the dual of H^(f2) via 



4, V0GH^(i7), 

J n Jdn 



(3.4) 



where (• | •) denotes the duality action on H^(i7) from its dual. We call D^w the 
generalised Hessian of w, and assuming that the coefficient tensor A is in C'^(]7)'*^'', 
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for the product with a distribution to make sense, we now seek u £ Hj!j(f2) such 
that Vu|g^2 ^ H~^/^(I2) and whose generaUsed Hessian satisfies 

(A:D2w|0) = (/,0) V(/)eHi(l2). (3.5) 
3.2. Finite element discretisation and finite element Hessian. We discretise 



(3.5) for simplicity with a standard piecewise polynomial approximation for test 
and trial spaces for both problem variable, and auxiliary (mixed-type) variable, 
H[U]. Let ^ be a conforming, shape regular triangulation of SI, namely, ^ is a 
finite family of sets such that 

(1) K £ implies K is an open simplex (segment for d = 1, triangle for d — 2, 
tetrahedron for d — 3), 

(2) for any K, J £ 3^ we have that X n J is a full subsimplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) [jK^^K^n. 

We use the convention where ft, : — > IR denotes the meshsize function of , i.e., 

h{x) : 

We introduce the finite element spaces 

\ := {<S> £ ll\n) : $|a' e P^'V/f e and $ e C"(/2)} , (3.7) 

V:=VnHj(J7), (3.8) 

where P*^ denotes the linear space of polynomials in d variables of degree no higher 
than a positive integer k. We consider p > 1 to be fixed and denote by N := dimV 
and N dimV. The discretisation of problem then reads: Find {U,H[U]) £ 

(iJ[C/],$) = - / VC/® V$+ / VC/®n<I> V$eV, 

Jn Jan (3.9) 

(A:i?[C/],*) = (/,«') V*eV. 



h{x) := m.ayihK- (3.6) 

K3x 



For an algebraic formulation of (3.9) we refer the reader to Lakkis and Pryer 2011 
§2]. Note that this discretisation can be interpreted as a mixed method whereby the 
first (matrix) equation defines the finite element Hessian and the second (scalar) 



equation approximates the original PDE (3.2) 



3.3. Two discretisation stategies of ( 1.1 ). The finite element Hessian allows us 



two discretisation strategies. The first strategy, detailed i n i|4} consists in applying 



Newton first to set-up (2.13) and then using the NVFEM (3.9) to solve each step. A 
second strategy becomes possible, upon noting that given U £ V the finite element 
Hessian H[U] is a regular function]^ which the generalised Hessian D^U might fail 
to be. This allows to apply nonlinear functions such as F to H[U] and consider 
the following fully nonlinear finite element method (FNFEM) 

(i?[C/], V[/«)V$+ / V[/®n$ V$eV, 

Jn Jan (3.10) 

(i^(J?[C/]),*) = (/,*) V^'eV- 

Of course, in order to solve the second equation, a finite-dimensional Newton 
method may be necessary (but this strategy leaves the door open for other nonlin- 
ear solvers, e.g., fixed point iterations). A finite element code based on this idea 
will be tested in ^ to solve the Pucci equation. 



generalised function ti is a regular function, or just regular, if it can be represented by a 
Lebesgue measurable function / £ 1}°'^ such that {v\<j>) = fif) for all £ CQ°(i7). We follow 
the customary and harmless abuse in identifying v with /. 
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In summary the finite element Hessian allows both paths in the following dia- 
gram: 

Newton 



fully nonlinear PDE (l.ll 



nonvariational linear PDE's (2.131 



fully nonlinear FE discretization (3.101 



Newton 



discrete linear 1 
discrete linear 2 



(3.11) 



Although the diagram in (3.11) does not generally commute, if the function F is 



algebraically accessible, then it is commutative. By "algebraically accessible" we 
mean a function that can be computed in a finite number of algebraic operations or 
inverses thereof. In this paper, we use only algebraically accessible nonlinearities, 
but, in principle assuming derivatives are available, our methods could be extended 
to algebraically inaccessible nonlinearities, such as Bellman's (or Isaacs's) operators 
involving optimums over infinite families, e.g.. 



inf Lc 



where {Lq 



F{M) 
a e £/} (or {La,i3 



M (orF(M):= inf supi„^:Af), (3.12) 

: (a, /?) e ^ X ^}) is a family of elliptic operators. 



4. The discretisation of unconstrained fully nonlinear PDEs 

In this section we detail the application of the method reviewed in f|3] to the 
fully nonlinear model problem (l.ll. Many fully nonlinear elliptic PDEs must be 



constrained in order to admit a unique solution. For example the Monge-Ampere- 
Dirichlet is elliptic and admits a unique solution in the cone of convex (or concave) 
functions when / > (or / < 0, respectively). Before we turn our attention to the 
more complicated constrained PDE's in S|5]and we illustrate the Newton-NVFEM 
method in the simplest light. In this section we study fully nonlinear PDEs which 
have no such constraint. 

4.1. Assumption (unconditionally elliptic linearisation). We assume, in this 
section, that the Newton-step linearisation (2.13) is elli ptic. For this assumption to 
hold, it is sufficient to assume uniform ellipticity, i.e., with = Sym(IR''^'^) 
and /i > independent of M . 

4.2. The Newton-NVFEM method. Suppose we are given a BVP of the form, 
finding u G H^(i7) n Ho(/?) such that 



jy[u] ^ FiB^u) - / = 



in i7. 



(4.1) 



which satisfies Assumption |4.1 



Upon applying Newton's method to approximate the solution of problem (4.1) 



we obtain a sequence of functions (M")ngiMo solving the following linear equations 
in nonvariational form, 

7V(D2m"):D2u"+i =5(d2u") 



where 



N{X) :=F'(X), 
g{X) := f-F{X) + F'iXy.X. 



(4.2) 

(4.3) 
(4.4) 
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The nonhnear finite element method to approximate (4.2) is: given an initial 



guess [/" := now° for each n € Nq find (C/"+i, JJ[[/"+i]) e V x V''^'' such that 

(iJ[C/"+^],$) + / VC/"+^ (g) V$- / VC/"+^(g)n<I> = V$eV 

Jf2 Jan (4.5) 

and (7V(i?[[/"]):Jf[J7"+i], {g{H[U"]), *) V * G V. 



4.4. Numerical experiments: a simple example. In this section we detail 



numerical experiments aimed at demonstrating the application of (4.5) to a simple 
model problem. 



4.5. Example (a simple fully nonlinear PDE). The first example we consider 
is a fully nonlinear PDE with a very smooth nonlinearity. The problem is 



^[u] := sin (Am) + 2Au - / = in i7, 

u = on dSl. 

which is specifically constructed to be uniformly elliptic. Indeed 

F\D^u) = (cos (Au) + 2) I . 



(4.6) 



(4.7) 



which is uniformly positive definite. The Newton linearisation of the problem is 
then: Given u^, for n € Nq find such that 



(cos (Alt") + 2) I: D^(m"+i - u") = / - sin (Au") - 2Aw" 
and our approximation scheme is nothing but |4.5| with 



N{X) = (cos (trace X) + 2)I 
g{X) = / — sin (trace X) — 2 trace X. 



(4.8) 



(4.9) 
(4.10) 



Figure [T] details a numerical experiment on this problem when d — 2 and when 
il = [—1, 1]^ is a square which is triangulated using a criss-cross mesh. 



4.6. Remark (simplification of Example 4.5). Example 4.5 can be simplified 
considerably by noticing that 



/ trff[L/]$= / (V[/)'^V$ V$eV. 



(4.11) 



This coincides with the definition of the discrete Laplacian and makes the NVFEM 
coincide with the standard conforming FEM. This observation applies to all fully 



nonlinear equations, with nonlinearity of the form (1.1) with 



F{M) := a(trM), 



(4.12) 



for some given a. This class of problems, can be solved using a variational finite 
element method and can be used for comparison with our method. Note that in 



Jensen and Smears 2012 the authors use this together with a localisation argu- 



ment in order to prove convergence of a finite element method for a specific class 
of Hamilton-Jacobi-Bellman equation. Their method coincides with ours, for an 
appropriate choice of quadrature. 
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Figure 1 

ing / appropriately such that u{x) 



Numerical experiments for Example 4.5 

( 



Choos- 
We use 



exp^— 10 \ x 

an initial guess = and run the iterative procedure until 
- J7"|| < 10-^ setting U := the final Newton iterate 
of the sequence. Here we are plotting log-log error plots together 
with experimental convergence rates for L2(r2),H\i7) error func- 
tionals for the problem variable, U ^ and an L2(J7) error functional 
for the auxiliary variable, H\U]. Notice that there is a "supercon- 
vergence" of the auxiliary variable for both approximations. 
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(a) Taking V to be the space of piecewise lin- (b) Taking V to be the space of piecewise qua- 

ear functions on (p = 1). Notice that dratic functions on i7 (p = 2). Notice that 

||«-i7*^|| = 0(h2), |„_[/A/|^ = o(fc) and = 0{h^), |m-C/"|j = 0(h?) and 

D2«-//[C/A-f]|| =0{h°-^). \\-D'^u- H[U^']\\ =0(hi-5) 



4.7. Example (nonvariational example). This is a simple example where the 
variational trick mentioned in Remark |4.6| cannot be applied. We fix d = 2 and 
consider the problem 



J/[u] := {diiuf' + {d22uf + diiu + d22U - / = in i7 

It = on dfl. 



(4.13) 



The approximation scheme is then (4.5) with 

1 



N{X) := 



iXl^ + 

g[X) :=/ + 2(X?i+Xy . 

Figure [2] details a numerical experiment on this problem in the case d = 2 and 
J? = [— 1, 1]^ triangulated with a criss-cross mesh. A similar example is also studied 



(4.14) 
(4.15) 



Davydov and Saeed 2012 Ex 5.2] using Bohmers method. 



5. The Monge-Ampere-Dirichlet problem 

In this section we propose a numerical method for the Monge-Ampere-Dirichlet 
(MAD) problem 

det D^M = f in n 

u ^ g on dfi. 

Our numerical experiments exhibit robustness of our method when computing 
(smooth) classical solutions of the MAD equation. Most importantly we noted 
the following facts: 
(i) the use of elements with p > 2 is essential as do not work, 



(5.1) 
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Figure 2 

ing / appropriately such that u{x) 



Numerical experiments for Example 4.7 

( 



Choos- 
We use 



exp^— 10 \ x 

an initial guess u'^ = and run the iterative procedure until 
- J7"|| < 10-^ setting U := the final Newton iterate 
of the sequence. Here we are plotting log-log error plots together 
with experimental convergence rates for L2(r2),H\i7) error func- 
tionals for the problem variable, U ^ and an L2(J7) error functional 
for the auxiliary variable, H\U]. Notice that there is a "supercon- 
vergence" of the auxiliary variable for both approximations. 
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(a) Taking V to be the space of piecewise lin- (b) Taking V to be the space of piecewise qua- 
ear functions on (p = 1). Notice that dratic functions on (p = 2). Notice that 

0{h) and ||u-(7^^|| = 0{h?), |m-C/"|j = 0(h?) and 



\u-U'"\\ = 0(h2), \u-U'^'\^ 



iB'^u- H[U'^'']\\ =0(hi-3) 



(ii) the convexity of the Newton iterates is conserved throughout the computa- 



tion, in a similar way to the observations in Loeper and Rapetti 2005b 



where the authors prove this convexity-conservation property. 
Our observations are purely empirical from computations, which leaves an interest- 
ing open problem of proving this property. 



5.1. Remark (the MAD problem fails to satisfy Assumption 4.1). To clar- 
ify Assumption 4.1 for the MAD problem (5.1 1, in view of the characteristic expan- 



0(e2 



(5.2) 



sion of determinant if X, 1" e Sym (R''^'^) 

det(X + eY) = dot X + eCof (X):Y 

where Cof X is the matrix of cofactors of X. Hence 

F'(X) = CofX. (5.3) 

This implies that the linearisation of MAD is only well posed if we restrict the class 
of functions we consider to those u that satisfy 



I^Cof D2m^> A|||^ V^e 



(5.4) 



for some (it-dependent) A > 0. Note that (5.4) is equivalent to the following two 
conditions as well 

(5.5) 



u is strictly convex. 

Loeper and Rapetti| (2005a have shown that for the continuous (infinite dimen- 

given an strictly convex initial guess 



2.4 



sional) Newton method described in 
each iterate m" will be convex. It is crucial that this property is preserved at the 
discrete level, as it guarantees the solvability of each iteration in the discretised 
Newton method. For this it the right notion of convexity turns out to be the finite 
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element convexity as developed in Aguilera and Morin 2009 . In Pryer 2010 , an 



intricate method based on semidefinite programming provided a way to constrain 
the solution in the case of elements. Here we observe that the finite element 
convexity is automatically preserved, provided we use or higher conforming ele- 
ments. 



5.2. Newton's method applied to Monge Ampere. In view of (5.3) it is clear 
that 

D^[u]v = Cof D^u: D\. (5.6) 
Applying the methodology set out in i|4]we set 

Ar(D2u") = Cof D^u", (5.7) 
g{ D^ti") = / - det D^u" + Cof D^u": D^u", (5.8) 

5.3. Remark (relating cofactors to determinants). For a generic (twice dif- 
ferentiable) function v it holds that 



ddet B^v = Cof B^v.B^ 



(5.9) 



Using this formulation we could construct a simple fixed point method for the 
Monge- Ampere equation. In view of Remark |5.3| ff can be further simplified 

g( D^u") = / - det D^u" + Cof D^u": D^m" 

(5.10) 

= /+ (d- l)det D^u". 

Newton's method reads: Given u'^ for each n e No find such that 

7V(D2u"):D2m"+i = .g(D2w"). (5.11) 

5.4. Numerical experiments. In this section we study the numerical behaviour 
of the scheme presented in Definition |4.2| applied to the MAD problem. 

We present a set of benchmark problems constructed from the problem data 
such that the solution to the Monge-Ampere equation is known. We fix fi to be 
the square S — [—1,1]^ or [0,1]^ (specified in the problem) and test convergence 
rates of the discrete solution to the exact solution. 

Figures [3]-[6] details the various experiments and shows numerical convergence 
results for each of the problems studied as well as solution plots, it is worthy of 
note that each of the solutions seems to be convex, however this is not necessarily 
the case. They are all though finite element convex Aguilera and Morin 2009|. In 



each of these cases the Dirichlet boundary values are not zero. The implementation 
of nontrivial boundary conditions is described in Lakkis and Pryer 2011 §3.6] or 
m more detail in [Pryer[ |2010l §4.4]. 



5.5. Remark (choosing the "right" initial guess). As with any Newton method 
we require a starting guess, not just for U° but also of H[U'^]. Due to the mild non- 
linearity with the previous example an initial guess of f/" = and H[U'^] = was 
sufficient. The initial guess to the MAD problem must be more carefully sought. 

Since we restrict our solution to the space of convex functions, it is prudent 
for the initial guess to also be convex. Moreover we must rule out constant and 
linear functions over i7, since the Hessian of these objects would be identically 
zero, destroying ellipticity on the initial Newton step. Hence we specify that the 



initial guess to (5.11 ) must be strictly convex. Rather than postprocessing the finite 



element Hessian from a initial project (although this is an option) to initialise the 
algorithm we solve a linear problem using the nonvariational finite element method. 
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Following a trick, described in Dean and Glowinski 2003 
standard V-finite element approximation of such that 

Au° = 2V7 in ^ 
u'^ — g on dfi. 



, we chose C/° to be the 

(5.12) 
(5.13) 



5.6. Remark (degree of the FE space). In the previous example the lowest 
order convergent scheme was found by taking V to be the space of piecewise linear 
functions {p = 1). For the MAD problem we require a higher approximation power, 
hence we take V to be the space of piecewise quadratic functions, i.e., p — 2. 

Although the choice of p = 1 gives a stable scheme, convergence is not achieved. 
This can be characterised by Aguilera and Morin 2009 Thm 3.6] that roughly 



says you require more approximation power than what piecewise linear functions 
provide to be able to approximate all convex functions. Compare with Figure [5j 



Figure 3. Numerical results for the MAD problem on the square 
5* = [— 1, 1]^. We choose the problem data / and g appropriately 
such that the solution is the radially symmetric function u{x) = 

exp^|x|^/2y We plot the finite element solution together with a 

log-log error plot for various error functionals as in Figure [ij Note 
for p = 2 the L2(i7) error rate of convergence is suboptimal, this 
is in agreement with the numerical examples produced in |Brenner| 



et al. 2011b 





(a) The FE approximation to the function (b) Log-log error plot for Lagrange FEs. 
u{x) = exp( J^). 



5.7. Nonclassical solutions. The numerical examples given in Figures |3}|6] both 
describe the numerical approximation of classical solutions to the MAD problem. 
In the case of Figure [sju e C°°(7?) whereas in Figure [eju G C°°{fl) n C°(Q). We 
now take a moment to study less regular solutions, i.e., viscocity solutions which 
are not classical. In this test we the solution 

u{x) = \x\^" (5.14) 



FEM FOR NONLINEAR ELLIPTIC PROBLEMS 



13 



Figure 4. A FE-convexity test for the numerical example given 
injsj We plot det H[U] together with the principle minor of H[U]. 




J 



(a) The principal minor of H[U], an approxi- 
mation to the Hessian of the function u{x) = 
expfWiy 




(b) The determinant of H[U]. 



Figure 5 . Numerical results for the MAD problem on the square 
S = [—1, 1]^. We choose the problem data / and g appropriately 
such that the solution is the radially symmetric function u(x) = 



exp|^|a;|'^/2 




(a) The FE approximation to the function (b) The error M—f7 plotted as a function over i?. 

Note the FE approximation does not converge 
• ' in this case, see Remark|5.6| 



uix) = expl 



for a e (1/2,3/4). The solution u{x) ^ W{[2). In Figures [7||9| we vary the 
value of a and study the convergence properties of the method. We note that 
the method fails to find a solution for a < 1/2. Finally in Figure 10 we conduct 
an adaptive experiment based on a gradient recovery aposteriori estimator. The 
recovery estimator we make use of is the Zienkiewicz-Zhu patch recovery technique 
see ?, 



Fryer 2010 §2.4] or Ainsworth and Oden 2000[ §4] for further details. 



14 OMAR LAKKIS AND TRISTAN FRYER 



Figure 6. Numerical results for the MAD problem on the square 
£'=[—1,1]^. Choosing / and g appropriately such that the solu- 

1 /2 

tion is u{x) = — (2 — — X2) ■ Note the function has singular 
derivatives on the corners of S. We plot the finite element solution 
together with a log-log error plot for various error functionals as 
in Figure [1] 



(a) The FE approximation to the function (b) Log-log error plot for Lagrange FEs. 
u{x) = -{2-xl-x1)^^^. 

Figure 7. Numerical results for the MAD problem on the square 
S = [—1,1]^. Choosing / and g appropriately such that the solu- 
tion is u{x) = \x\ with a = 0.55. Note the function is singular 
at the origin. We plot the finite element solution together with a 
log-log error plot for various error functionals as in Figure [l] 
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Figure 8. Numerical results for the MAD problem on the square 
= [— 1, 1]^. Choosing / and g appropriately such that the solu- 
tion is u{x) = |a;|^", with a = 0.6. Note the function is singular 
at the origin. We plot the finite element solution together with a 
log-log error plot for various error functionals as in Figure [l] 




Figure 9. Numerical results for the MAD problem on the square 
S = [—1,1]^. Choosing / and g appropriately such that the solu- 
tion is u{x) — |a;|^", with a = 0.7. Note the function is singular 
at the origin. We plot the finite element solution together with a 
log-log error plot for various error functionals as in Figure [T] 
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(a) The FE approximation to the function (b) Log-log error plot for P-^ Lagrange FEs. 
u{x) = Ixl^", with a = 0.7. 
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Figure 10. Numerical results for a solution to the Monge- 
Ampere-Dirichlet equation with / and g appropriately such that 
the solution is u{x) = with a — 0.55. We choose p = 2, 

and use an adaptive scheme based on Z-Z gradient recovery. The 
mesh is refined correctly about the origin. Note that when dim V = 
20, 420 the adaptive solution achieves \\u— C/*^ || « 0.0078, the uni- 
form solution given in Figured satisfies \\u ~ ?7*^ || ~ 0.18 using the 
same number of degrees of freedom. Using the adaptive strategy 
both II u — ?7*^|| and \u — W^'l^ converge like 0(iV^^). 




(a) Adaptive mesh 



(b) Log-log error plot for Lagrange FEs. 
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6. PUCCl'S EQUATION 

In this section we look to discretise the nonlinear problem, in this case Pucci's 
equation as a system of nonlinear equations. Pucci's equation arises as a linear 
combination of Pucci's extremal operators. It can nevertheless be written in an 
algebraically accessible form, without the need to compute the eigenvalues. 



6.1. Definition (Pucci's extremal operators Caffarelli and Cabre 1995| ). 

Let TV G Sym [R'^^'^) and (t{N) be it's spectrum, then the extremal operators are 

^(AT) := ^ aA. = 0, (6.1) 

with ai e R. The maximal (minimal) operator, commonly denoted (^~), 
has coefficients that satisfy 

< ai < • • • < a„ (ai > • • • > a„ > 0) (6.2) 

respectively. 



6.2. The planar case and uniform ellipticity. In the case d = 2 the normalised 
Pucci's equation reduces to finding u such that 



aX2 + Ai = 



(6.3) 



where N :— D^u. Note that if a = 1 (6.3 1 reduces to the Poisson-Dirichlet 
problem. This can be easily seen when reformulating the problem as a second order 
PDE Dean and Glowinski 2005 . Making use of the characteristic polynomial, we 
see 



A,; = 



Au±({Auf - Adet D'^uj 



1/2 



z = 1,2. 



Thus Pucci's equation can be written as 

= (a + 1) Au + (a - 1) ({Auf - 4 det D^w 



1/2 



(6.4) 



(6.5) 



which is a nonlinear combination of Monge- Ampere and Poisson problems. However 
owing to the Laplacian terms, and unlike the Monge-Ampere-Dirichlet problem, 
Pucci's equation is (unconditionally) uniformly elliptic for 



(trX) -4detX > VXe 



(6.6) 



The discrete problem we use is a direct approximation of (6.5 ), we seek(C/, H[U]) 
such that 



^ (^{a + 1) tr H[U] +{a-l) ((tr H[U]f - 4 dot H[U]j ^'"^ j $ 



= 



(ff[[/],*) = - / V* 



dn 



VL/^n* V($,*)eVxV. 



(6.7) 
(6.8) 



The result is a nonlinear system of equations which was solved using a algebraic 
Newton method. 
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Figure 11. Numerical results for a classical solution to Pucci's 



equation (6.9). As with the case of the MAD problem we choose 
p = 2. We use a Newton method to solve the algebraic system until 



the residual of the problem (see Kelley, 1995 c.f.]) is less than 
10""'^° (which is overkill to minimise Newton error effects). We 
plot log-log error plots with experimental orders of convergence, 
for various norms and values of a. 
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6.3. Numerical experiments. We conduct numerical experiments to be com- 
pared with those of Dean and Glowinski 2005 . The first problem we consider is a 



classical solution of Pucci's equation (6.3). Let x ={x,yy, then the function 

,(l-a)/2^ 



uix) = -((ix + if+iy + if) 



(6.9) 



solves Pucci's equation almost everywhere away from {x,y) = (—1,-1) with g 
u\gQ. Let 3^ be an irregular triangulation of J7 = [—0.95,1]^. In Figure 
detail a numerical experiement considering the case a G [2, 5]. 
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we 



We also conduct a numerical experiment to be compared with Oberman 2008 
In this problem we consider a solution of Pucci's equation with a piecewise defined 
boundary. Let H ~ [—1,1]^ and the boundary data be given as 



when \x\ > | and \y\ 
otherwise. 



> i 



(6.10) 
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Figure 12. Numerical results for a solution to Pucci's equation 
with a piecewise defined boundary condition (6.10). We choose 
p = 2, use a Newton method to solve the algebraic system until the 
residual of the problem is less than 10^^°. We plot the solution for 
various values of a as well as a cross section through the coordinate 
axis. Notice that the solution becomes extremely badly behaved 
as a increases. 





Figure [T2| details the numerical experiment on this problem with various values of 

Q. 



Since the solution to the Pucci's equation with piecewise boundary (6.10) is 
clearly singular near the discontinuities we have also conducted an adaptive exper- 
iment based on a gradient recovery aposteriori estimator (as in ^5.7). As can be 
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Figure 13. Numerical results for a solution to Pucci's equation 
with a piecewise defined boundary condition (6.10). We choose 



p = 2, and use an adaptive scheme based on Z-Z gradient recovery. 
The mesh is refined correctly about the jumps on the boundary. 




j 



(a) Finite element solution for a = 3 



(b) Adaptive mesh 



seen from Figure [13] we regain qualitively similar results using far fewer degrees of 
freedom. 



7. Conclusions and outlook 

In this work we have proposed a novel numerical scheme for fully nonlinear and 
generic quasilinear PDEs. The scheme was based on a previous work for nonvaria- 



tional PDEs (those given in nondivergence form) Lakkis and Pryer 2011 



We have illustrated the application of the method for a simple, non physically 
motivated example, moving on to more interesting problems, that of the Monge- 
Ampere equation and Pucci's equation. 

For Pucci's equation we numerically showed convergence and conducted experi- 
ments which may be compared with previous numerical studies. 

We demonstrated for classical solutions to the Monge-Ampere equation the 
method is robust again showing numerical convergence. For less regular viscoc- 
ity solutions we have found that the method must be augmented with a penalty 



term in a similar light to Brenner et al. 2011b 



We postulate that the method is better suited to a discontinuous Galerkin frame- 
work which is the subject of ongoing research. 
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